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January  10,  1967 


The  FORTRAN  IV  computer  cudc  lor  the  UNIV AC  1108  described 
herein  is  as  it  existed  on  12/1/66.  The  code  has  been  in 
continuous  development  for  1  year  and  in  its  presented  form 
has  been  applied  successfully  by  General  Atomic  to  the  kind 
of  problems  discussed  later  in  this  report.  However,  the 
development  and  improvement  of  the  code  are  being  continued, 
so  that  duplication  of  results  (or  even  close  agreement) 
between  problems  run  with  the  code  as  published  and  the  code 
as  it  existed  either  before  or  after  this  time  is  not  necessarily 
to  be  expected. 

General  Atomic  has  exercised  due  care  in  preparation,  but 
does  not  warrant  the  merchantability,  accuracy,  and  complete¬ 
ness  of  the  code  or  of  its  description  contained  herein.  The 
complexity  of  this  kind  of  program  precludes  any  guarantee  to 
that  effect.  Therefore,  any  user  must  make  his  own  determina¬ 
tion  of  the  suitability  of  the  code  for  any  specific  use,  and  of 
the  validity  of  the  information  produced  by  use  of  the  code. 
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ABSTRACT 


A  program  is  presented  to  compute  fallout  fission  product  absorption 
in  an  expanding,  cooling,  uniform  field  of  gas  and  fallout  particles  where 
the  rate  of  fission  product  absorption  is  controlled  by  fission  product 
surface  concentrations  as  given  by  Henry's  law  constants  and  diffusion  of 
these ussion  products  into  the  fallout  particles.  The  calculations  are  made 
a  nuclide  cnainat  a  time,  employing  nuclear  device  and  fission  product 
parameters.  Program  output  includes  average  concentration  versus 
particle  size  for  each  absorbed  fission  product  at  a  preselected  lower 
temperature  and  the  amount  of  that  fission  product  in  the  gas  phase,  with 
an  option  for  the  calculation  of  fission  product  radial  profiles  versus 
particle  size.  As  a  demonstration  of  this  model,  an  initial  computer  study 
of  fallout  formation  phenomenology  is  presented. 
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SUMMARY 


The  computer  code  described  in  this  report  represents  a  model  of 
fallout  fission  product  absorption,  the  fallout  particles  being  distributed 
uniformly  in  a  cooling,  expanding  field  of  gas.  Rates  of  fission  product 
absorption  are  governed  by  fission  product  surface  concentrations  as 
given  by  Henry's  law  constants  and  by  diffusion  into  the  fallout  particles. 

It  is  assumed  that  the  particles  neither  agglomerate  nor  convect  fission 
products  during  the  time  stepped  calculational  portion  of  the  program. 

The  absorption  of  fission  products  is  considered  for  one  decay  chain  at 
a  time  using  detonation  yields,  initial  fission  product  yeilds,  half-lives, 
Henry's  law  constants,  concentration  independent  diffusivities,  and  a 
fallout  particle  size  versus  quantity  distribution.  Output  consists  of 
surface  concentration  versus  temperature  for  each  absorbed  fission 
product,  average  concentration  versus  particle  size  at  a  final  lower 
temperature,  and  the  amount  of  that  fission  product  remaining  in  the  gas 
phase.  Optionally,  radial  profiles  of  fission  product  concentration  versus 
particle  size  at  this  temperature  are  calculated. 

Fission  product  distribution  for  mass  chains  95  and  137  in  fallout 
using  this  code  are  described.  One  of  the  important  findings  derived  from 
these  computer  studies  is  an  understanding  of  how  this  model  considers 
fractionation  to  vary  with  fallout  particle  size  distribution. 


v 


1.  INTRODUCTION 


Since  the  phenomenological  discovery  of  fractionation  of  fission 
product  elements  in  recovered  fallout,  the  importance  of  the  chemical 
properties  of  the  fission  products  and  fallout  particles  to  the  description 
of  fallout  formation  has  been  recognized.  In  his  treatise  on  fallout, 

Miller^  has  attempted  to  synthesize  fractionation  effects  by  allowing 

o 

fission  products  to  equilibrate  with  fallout  until  they  have  cooled  to  1400  C 
and,  thereafter,  to  surface  condense  on  the  fallout.  This  scheme  has 
been  considered  an  oversimplification  and  phenomenologically  unsatisfactory. 
General  Atomic  has  taken  the  position  that  condensed  state  diffusion  of 
fission  products  in  glassy  silicate  fallout  particles  can  be  considered  to 
control  the  rate  of  fission  product  absorption.  This  report  is  a  description 
of  a  calculational  program  written  in  FORTRAN  IV  for  the  UNIVAC  1108 
in  which  the  surface  of  fallout  is  considered  to  be  in  equilibrium  with  the 
gas  phase,  but  bulk  fallout  phase  absorption  of  fission  products  is  diffusion 
controlled.  In  addition,  fallout  formation  information  derived  from  some 
computer  experiments  employing  this  model  is  presented  in  Appendix  B. 
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2.  GENERAL  DESCRIPTION 


The  problem  consist*  of  determining  concentrations  of  radioactive 
isotopes  in  fallout  particles.  Absorption  of  fission  products  is  considered 
a  decay  chain  at  a  time  where  a  chain  consists  of  a  number  of  isotopes 
generally  of  the  same  mass.  The  last  isotope  in  the  chain  is  assumed  to 
be  stable  for  the  time  periods  being  considered.  The  fallout  particles  are 
assumed  to  be  non-agglomerating,  internally  non-convecting,  refractory, 
spherical  glassy  silicate  bodies  of  various  sizes  uniformly  distributed 
throughout  a  uniform  temperature-concentration  cloud  of  particles,  air, 
and  fission  products.  Immediately  after  the  detonation  takes  place,  tem¬ 
peratures  are  very  high  and  chemistry  in  the  cloud  is  unimportant. 
However,  at  a  time  when  the  cloud  has  cooled  sufficiently  to  consider 
chemistry,  this  program  can  be  applied.  From  this  point,  as  time  pro¬ 
gresses  (and  the  temperature  of  the  cloud  further  decreases)  the  isotopes 
are  assumed  to  surface  condense  according  to  Henry's  law  and  subsequently 
to  be  diffused  into  the  particles.  Concurrently,  the  isotopes  decay 
from  one  to  another  at  rates;  dependent  upon  their  respective  half-lives. 

At  a  sufficiently  low  temperature,  the  diffusion  process  is  considered  to 
terminate  leaving  the  isotopes  distributed  among  and  within  the  particles 
in  some  manner.  The  purpose  of  the  program  is  to  describe  that  dis¬ 
tribution  in  terms  of  surface  co.  centrations,  average  concentrations,  and 
concentrations  throughout  the  particles  (profiles).  To  find  these  various 
concentrations,  it  suffices  to  calculate  gas  and  condensed  phase  mass 
balances  in  which  surface  concentrations  of  each  isotope  are  determined 
as  a  function  of  time. 

Time -temperature  stepping  is  employed  in  the  program,  and  mass 
balances  are  determined  at  each  time  step.  In  order  to  determine  the  mass 


3 


balances,  it  is  necessary  to  use  input  nuclear  cloud  parameters  and 

chemical  information.  Time  scaling  functions  for  temperature  and  for 

cloud  volume  and  a  fallout  mass  representation  are  employed  according 

to  Miller^  ^  but  can  be  changed  to  fit  other  circumstances.  Fission  product 

yields  can  be  taken  from  Crocker,  *2)  and  fission  product  decay  rates  can 

(3) 

be  taken  from  Crocker,  Scheldt,  and  Connors.  Henryk  law  constants 

(4) 

for  this  program  have  been  presented  by  Norman.  A  similar  report 
on  diffusion  coefficients  for  this  program  is  planned. 

Surface  concentrations  are  obtained  directly  from  the  mass  balance 
calculations.  Other  output  parameters  are  calculated  from  the  surface 
concentration  values  and  the  degree  of  diffusion  that  has  been  calculated. 

The  computation  of  concentration  profiles  of  the  various  sired  particles 
at  a  final  temperature  is  optional.  Another  option  allows  calculation  of 
late  distributions  at  such  a  time  that  all  remaining  gaseous  short-lived 
isotopes  have  decayed  and  condensed.  In  this  option,  it  is  assumed  that 
the  isotopes  are  deposited  on  the  particles  according  to  surface  area. 

Surface  concentrations  are  calculated  at  a  finite  number  of  points 
assuming  that  an  actual  time -temperature  distribution  can  be  closely 
approximated  by  a  step  function  as  shown  in  Fig.  1.  The  exact  nature  of 
this  step  function  depends  on  the  starting,  final,  and  incremental  tem¬ 
peratures  (seeSection  4. 1,  item  4).  Thus,  all  temperature  (and  time) 
dependent  quantities  are  assumed  to  be  step  functions  as  well;  in  particular, 
volumes  and  diffusion  constants  and  Henry’s  law  constants  for  each  isotope 
remain  constant  for  each  temperature  interval.  Fission  product  surface 
concentrations  are  also  assumed  constant  during  any  time  interval. 
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3.  METHOD 


To  determine  the  surface  concentration  of  a  given  isotope  at  some 
temperature,  the  surface  concentrations  of  that  isotope  and  all  its  pre¬ 
cursors  at  all  preceding  temperatures  must  be  employed,  since  some 
fraction  of  each  precursor  will  decay  into  the  next  isotope  an*  previous 
surface  concentrations  will  affect  the  total  surface  concentration  at  the 
specified  temperature.  We  start  by  assuming  that  all  of  the  fission 
products  are  uniformly  distributed  throughout  the  cloud  at  some  sufficiently 
high  temperature.  To  find  surface  concentrations  at  the  first  temperature, 
the  following  expression  is  used: 


ciiviMjHn 

RT, 


“4 


D(#>w 
ill  I 


=  Y 


il 


(1) 


where  C.  ^  =  surface  concentration  (g/g)  of  isotope  i  at  time  t^  (sec) 

Vj  =  volume  of  nuclear  cloud  at  time  tj  (liters) 

M.  =  molecular  weight  of  isotope  i  (g/mole) 

H.  j  =  Henry's  law  constant  for  isotope  i  at  time  t  j  (atm/g/g) 

R  =  ideal  gas  constant  (liter  atm/°K  mole) 

Tj  =  temperature  at  time  t  j  (°K) 

Yfi  =  yield,  or  total  amount  of  isotope  i,  available  at  time  t^  (g) 

w  =  weight  of  soil  (g)  contained  in  particle  size  fraction  s  of 
radius  r8  (cm) 

d[®|  =  solution  of  the  diffusion  equation  as  if  diffusion  had  taken 

place  at  a  constant  surface  concentration  Cjj  (see  Section  4.  3) 
for  particle  of  s**1  radius,  with  the  6|]^tjterm  for  isotope  i 
assumed  constant  during  At}  (6i  1  Is  the  pertinent  diffusion 
coefficient (6  in  cm^/sec)).  The  diffusion  coefficients  are 
assumed  to  be  concentration  independent. 

^Suggested  self-consistent  units. 
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The  Yjj  value  is  computed  by  considering  the  decay  of  isotopes  and  the 
H^j  value  is  computed  from  characteristics  of  the  isotope  and  temperature 
(see  Section  4.  1 ). 

The  first  term  on  the  left  of  Eq  (1)  represents  the  amount  of  isotope 
i  in  the  gas  phase  at  time  1.  The  second  term  is  the  amount  in  the  solid, 
or  condensed,  phase.  For  each  s, 


(•) 

ill 


for  the  a1*1  particle  size,  where  is  the  average  concentration  during 

time  t^.  Thus 
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It  should  be  noted  that  in  this  program  an  approximation  system  is 
used  in  representing  the  decay  of  fission  products.  During  the  time 
intervals  used,  a  fission  product  is  assumed  to  decay  to  one  and  only  one 
product.  Thus,  decay  is  not  treated  as  a  dynamic  process  but  as  a  dis¬ 
crete  process  occurring  only  at  the  end  of  a  time  step  period.  A  nearly 
correct  fission-product  distribution  is  obtained  at  the  initiation  of  the 
program  by  stepping  the  time  in  many  small  increments  to  this  initial 
time  after  detonation.  The  time  increments  then  become  governed  by  the 
temperature  increments.  This  system  will  cause  some  deviation  from 
an  exact  decay  solution  but,  generally,  will  only  alter  the  final  answers 
in  a  negligible  fashion.  Of  course,  as  the  temperature  increment  approaches 
zero,  the  decay  system  will  approach  the  theoretical  limit. 
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To  find  surface  concentrations  at  the  second  temperature,  the 
expression  is 


Ci2V2MiH12 

RT, 


+  Cu  exp(-k1At1) 


i“'l 


Ui,  1,2s 


(1  '  ®xp(-ki-iAti,] 


w 

(i- 1 ).  1,  l;i.  2,2  s 


+  |Ci2  '  Cil  exp("kiAtl)mC(i-l)l  ^ 


exp(*k(i-l)At 


_w  =  Y 
2  s  i2 


(2) 


The  terms  on  the  left  side  of  Eq.  (2)  can  be  explained  as  follows: 

1.  The  first  term  is  due  to  the  amount  of  isotope  i  still  in  the  gas 
phase. 

2.  The  second  term  is  due  to  the  amount  of  i  surface  condensed 

associated  with  time  tj  decaying  for  Atj  and  being  diffused  for 

6. .At,  +  6.,  At,, 
xl  1  i2  2 

3.  The  third  term  is  due  to  the  amount  of  isotope  i- 1  (the  immediate 
precursor  at  tj  of  i)  surface  condensed  at  time  tj  decaying 

into  i  and  being  diffused  as  i-1  during  time  Atj  and  as  i  during 
time  At2> 

4.  The  fourth  term  represents  the  amount  of  i  being  diffused  within 
the  particles  as  the  result  of  the  surface  concentration  incre¬ 
ment  caused  by  the  temperature-time-volume  increment. 

The  D1 8  are  further  defined  in  Appendix  A. 

It  is  apparent  that  as  time  is  increased,  the  number  of  terms  in  the 
expression  for  C.^  becomes  quite  large.  In  fact,  it  can  be  shown 
(Appendix  A)  that  if  there  are  n-1  predecessor"  to  an  isotope,  the  number 
of  terms  required  to  find  C.^  is 
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where  0  <  r  <  i-n+1 
0  <  j  <  N-l 

The  decay  term  depends  on  the  position  of  the  isotope  i-r  in  the  chain 

and  on  the  time.  D.  and  D._.  are  diffusion  terms  which  also  depend 

in  i^) 

on  the  position  and  time.  (See  Appendix  A  for  details.  ) 

Once  surface  concentrations  have  been  found  for  all  times,  it 
remains  to  compute  average  concentrations  and  concentrations  at  various 
radii  within  the  particles.  Average  concentrations  can  be  obtained  by 
mathematically  combining  concentration  increments  according  to  an 
additivity  rule,  as  explained  in  Appendix  A.  Concentrations  at  any  point 
within  a  particle  can  be  found  in  an  analagous  way  using  relationships 
given  by  Crank.  ^ 


. 

m+r‘ 
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4.  SUBPROGRAM  DESCRIPTIONS 


This  section  describes  each  of  the  FORTRAN  IV  subroutines  used  by 
this  code.  Most  of  the  variables  used  are  stored  in  one  of  the  common 
blocks  described  in  Table  1. 

4.  1.  SUBROUTINE  INITAL 

This  routine  reads  the  data  cards  for  each  calculation  and  prepares 
the  input  data  for  the  main  program.  The  flow  is  in  approximately  the 
order  described  below. 

1.  Read  input  cards.  The  following  information  is  read  for  each 
calculation  with  the  internal  variable  name  indicated  in  each 
case.  (See  Appendix  C  for  card  formats.  ) 

a.  Length  of  the  chain,  LC  (<6) 

b.  Yield  of  device  used  for  calculating  fission  product  yields 
(in  kilotons  TNT),  YKT 

c.  Yield  of  device  used  for  calculating  the  time-temperature 
distribution  and  the  cloud  volume  (in  kilotons  TNT),  BKT 

d.  For  each  isotope  in  the  chain: 

(1)  Chemical  symbol,  NAME 

(2)  Molecular  weight  (in  g/mole),  WMOL 

4 

(3)  Initial  yield  (in  atoms/  1 0  fissions),  YI 

(4)  Half-life  (in  seconds),  HL 

(5)  Constants  used  in  calculating  Henry's  law  constants  and 
diffusion  coefficients  (see  item  8),  DC1,  DC2,  HC1,  and 
HC2 

e.  Starting  temperature,  final  temperature  and  incremental 

temperature  (in  °K),  HTMP,  ETMP,  and  TEMPV  (also 

referred  to  as  T8,  Tf,  and  AT) 
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TABLE  I 

VARIABLES  IN  COMMON  BLOCKS 


Variable 

Dimension! 

Description 

BLANK  Common:  Contained  in  NORMAN,  BCALC,  INITIAL 

YLD 

(6.40) 

Yield  of  each  ieotope  at  each  time  step 

NAME 

(6) 

Chemical  symbol  of  each  isotope 

WMOL 

(6) 

Molecular  weight  of  each  isotope 

DFC 

(6,40) 

Diffusion  coefficients  of  each  isotope 
at  each  time  step 

HEN 

(6,40) 

Henry's  law  constants  for  each  isotope  at 
each  time  step 

AK 

(6) 

Exponential  constant  associated  with 
each  half>life 

R 

(30) 

Radius  of  each  particle  size 

WT 

(30) 

Weight  of  soil  in  each  particle  size 

VOL 

(40) 

Volume  of  the  cloud  at  each  time  step 

T 

(40) 

Time  increments 

TEMP 

(40) 

Temperatures  at  each  time  step 

ATM 

(40) 

Times  at  each  step 

SI 

MG:  Contained  in  NORMAN,  FINAL,  BCALC,  INITIAL 

LC 

(i) 

Length  of  chain 

IR 

0) 

Number  of  particle  sizes 

NTOT 

(1) 

Number  of  time  steps 

NCS1 

(1) 

Case  identification 

NCS2 

(1) 

Case  identification 

TWT 

(1) 

Total  mass  of  soil 

ILC 

(1) 

Flag:  if  =  0,  do  not  do  FINAL  calculation 
if  i  0,  do  FINAL  calculation 

IDPTH 

(1) 

Number  of  radii  used  lor  profile  calculations 

Q 

(11) 

Relative  radii  for  profile  calculations 

Bl:  Contained  in  NORMAN,  BCALC,  DELTA 


DCY 

(6.  40, 40) 

Array  of  decay  constants  for  each  isotope  from 
each  time  step  to  each  other 

WN 

(6,  40) 

Array  of  concentration  ratios  times  weights  for 
each  isotope  at  each  time  step 

INLOOP:  Contained  in  NORMAN,  BCALC,  DELTA 

L 

(0 

Subscript  of  isotope  presently  under  consideration 
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f.  Number  of  particle  aizei,  IR  (<.30) 

g.  For  each  particle  size,  its  radius  (in  cm)  and  the  fraction  of 
soil  contained  in  particles  of  that  size,  R,  PW 

h.  The  total  mass  of  soil  (in  g/KT),  TWGT 

2.  Check  starting  temperature.  If  the  starting  temperature  is  so 
high  that  all  the  fission  products  are  uniformly  distributed 
throughout  each  particle,  that  temperature  is  lowered  by  one 
increment  and  the  procedure  is  repeated  until  a  temperature  is 
reached  where  this  is  not  the  case.  Calculations  are  initiated 
after  stepping  back  one  temperature  increment  or  to  the  starting 
temperature.  To  accomplish  this,  a  D  for  the  largest  particle 
radius  is  found  for  each  isotope  at  an  initial  time  step.  A  D  of  1 
implies  that  an  isotope  is  uniformly  distributed  throughout  all 
particles.  If  D  =  1  for  all  isotopes  and  particle  sizes,  the  process 
is  not  yet  diffusion  controlled;  thus,  if 

LC 

X)  Di  =  LC  ’ 

i=l 

all  D's  are  unity.  (Note  that  D<,1.  )  If  this  is  satisfied,  the 
temperature  is  reduced  and  the  procedure  is  repeated. 

3.  Calculate  w  ,  the  weight  of  soil  in  each  particle  size: 
wB  =  TWT  •  PWfl,  where  TWT  =  TWGT  -  BKT. 

(2) 

4.  Calculate  the  time-temperature  distribution.  Using  the  formula 

T  =  4.66  x  103  (BKT)’°*  010  exp(-0.546(BKT)"°'373t)  =  f(t)  ,  (6) 

times  and  incremental  times  are  calculated  to  correspond  to  each  specified 
temperature.  The  temperatures  are  calculated  from 


and  the  times  such  that 

Tn  =  f(tn>  * 
n  n 
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At  la  defined  by 
n  ' 

AV'‘1(v¥)-f'1(T„-T:)<*«F1«-  »  •  <7» 

5.  Convert  initial  yields,  Y^,  to  grams. 

Yt(g)  =  YI  •  WMOL  •  YKT  •  2.  3410261  x  10’5 

6.  Find  the  quantities  k^.  These  quantities  are  used  to  calculate 
decay  coefficients  where 

dlnm  *  +  At„+1  +  ' - +  Atm>l  m  *  n 


and 


=  exp(k  ,AO  m  =  n 
=  1  (a  computational  aid)  m  <  n 


ln(2) 

"  HLi  * 


7.  Calculate  yields  at  tg  where 


t  =  f  (T  )  . 

8  S 


Yields  at  t  are  approximated  by  the  following  scheme: 
8 


Let  Ah  =  min(HL)/10  when  min(HL)  is  the  minimum  half-life 

in  the  chain 


Let  t  =  n  (Ah),  and  let  i-1  specify  the  precursor  of 
isotope  i.  Then, 


Y.(Tl)=  Yi(0)exp(-k.Ah)  +  Yil(0)[l  -  exp(-ki_  jAh)]  (8) 

Yi*Tn>  =  YitTn.1)  ex?(-ki^h)  +  Yi.i(Tn.l^1  '  exp(_ki.lAh^  (9) 
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Equation  (9)  1>  repeated  for  ail  nand  i  until  n  ia  such  that 
Tn+1  >  S'  *'or  *aBt  BteP*  Ah  i*  changed  auch  that  Ah  = 
t  -  T^and  Eq.  (9)  ia  repeated.  If  the  original  Ah  >.  t^ ,  Ah 
ia  aet  equal  to  t  and  only  Eq.  (8)  ia  uaed. 

8.  Calculate  diffuaion  coefficienta  and  Henry'a  law  conatanta. 

(4) 

The  expressions 

log(H,  )  =  HC1.  -  HC2./T  (10) 

°  in  i  i  n 

and 

log(6in)  =  DClj  -  DC2i/Tn  (11) 

are  used  to  find  8  and  H  for  all  isotopes  and  temperatures. 

9.  Volumes  are  computed  for  all  temperatures  using  fhe 

.  (3) 

expression 

3 

V  =  |5.  69  x  102(BKT)1/3  exp(0.  17(BKT)"°' 3?3t  ]j  4“ 

10.  Input  and  computed  quantities  are  printed.  Input  ia  actually 
printed  as  it  is  read  and  calculated  quantities  as  they  are 
calculated.  (See  Appendix  D  for  sample  printout.  ) 

4.  2.  MAIN  ROUTINE  (NORMAN) 

The  main  routine  has  three  basic  functions: 

1.  Setup  and  initial  calculation  of  arrays  needed  (such  as  yields 
and  decay  coefficients)  for  the  main  loop. 

2.  Execution  of  the  main  loop  which  finds  surface  concentrations, 
average  concentrations,  and  profiles. 

3.  Writing  out  of  results. 

4.  2.  1.  Setup  and  Initial  Calculations 

The  main  routine  starts  by  calling  subroutine  INITAL,  which 

supplies  it  with  necessary  input  data.  Next  it  computes  the  array  of 

decay  coefficients,  d.  ,  as  defined  in  Section  4.  1,  item  6;  the  array  of 
7  inm  7 
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/ 


yield*  Y^.  where  Y^  is  the  yield  of  isotope  i  at  timen;  and  the  arrays  UN 
and  WN.  The  yield*  are  found  by  using  Eq.  (8).  Section  4. 1,  item  7, 
modified  as 


Y.  *  Y.  ,d.  ,  ,  +  Y.  ,  . (1  -  d.  .  .  ,)  . 

in  i»n-l  i.n-l»n-l  i- 1  # n- 1  i-l,n-l,n-l 


The  UN  and  WN  array*  are  defined  by 

UNin*=  E<6i-At-'r-) 


WN.  =  > 

ln  -fef 


r.  9  •  / 

in  n  s 


UN.  w 

ins  • 


where 


6.  =  diffusion  coefficient  of  Isotope  i  at  time  t 

in  ^  n 

At  =  n  time  interval 

n  th 

r  =  radius  of  •  particles 

•  th 

=  weight  of  soil  in  the  s  particle  size 

E(6.  At  ,  r  )  =  solution  of  diffusion  equation  for  sphere  of  radius  r  and 

in  A  S  m  m  «  «  1  _  8 

and  6,  At  term 
in  n 

Finally,  the  surface  concentrations  are  calculated  for  each  isotope  at 


time  tj  by  the  relation 


/  WMOL.H..V. 

I -  *  +  WN 

il  \  RT, 


a)  ■  yh 


4.  2.  2.  Main  Loop 

The  following  steps  are  taken  in  the  main  loop  for  each  isotope  and 
each  time: 

1.  Find  the  terms  due  to  decay  of  this  isotope  from  previous  times. 

2.  For  each  precursor  of  this  isotope  at  each  preceeding  time, 
find  all  possible  paths  for  decay,  and  for  each  such  path  find  a 
diffusion  term  (see  Section  4.  8)  and  a  decay  term  (see  Sectior  4.7), 

tUNln8  =  D\n; 

'See  Fig.  2 
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Fig.  2-*  Flow  chart,  main  loop  NORMAN  (Sheet  2  of  5) 
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For  each  desired  radios,  add  the  lurm 

PDELT  •  CNEp,  N22-  U  (DTERM.r.a>  -  HDTERM  I .  r.  •)) 

to  the  proper  com  entration  »om 


Fig.  2- -Flow  chart,  main  loop  NORMAN  (Sheet  3  of  5) 
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Fig.  2--Flow  chart,  main  loop  NORMAN  (Sheet  4  of  5) 
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Fig.  2- -Flow  chart,  main  loop  NORMAN  (Sheet  5  of  5) 
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forming  the  product  as  in  Eq.  (5)  and  adding  to  the  previous 

summed  contributions.  In  the  event  that  the  final  time  step 

has  been  reached,  add  the  appropriate  terms  to  the  average 

concentration  and  profile  point  sums. 

3.  Let  a  represent  the  total  of  the  sum  of  all  contributions  to  the 

average  concentration  expression  except  that  term  due  to 

absorption  diffusion  of  this  isotope  from  surface  concentration 

C.  .  Then  the  surface  concentration  is  found  from 
in 


C. 

in 


( 


WMOL.H.  V 
i  in  n 


RT 


n 


+  WN 


a  . 


4.  2.  3.  Output 


Certain  other  calculations,  as  follows,  are  made  while  printing 
out  results. 

1.  Calculation  of  the  amount  of  material  remaining  in  the  gas 
phase  for  each  isotope  at  each  time: 


GAS. 

in 


C.  H.  V  WMOL. 
in  m  n _ l 

RT 

n 


2.  Computation  of  the  amount  of  material  in  the  solid  phase 
(diffused),  per  isotope,  at  the  final  time: 

SOLID.  = 

where  A.g  is  the  average  concentration  of  isotope  i  in  particle 
size  s. 

3.  Computation  of  the  weighted  mean  average  concentration  for 
each  isotope  i: 

SOLID. 

A  - - - - 

i  (total  weight) 
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The  following  quantities  are  printed:* 

a.  T  ,  t  ,  and  C.  for  each  i,  n 

n  n  in 

b.  T  .  t  ,  and  GAS.  for  each  i.  n 

n  n  in 

c.  T„,  r  ,  and  A,  ,  N  being  the  final  time  and  for 

N*  is 

each  i,  s 

d.  A  for  each  i 

e.  C.  (r)  for  each  r,  i.  s,  where  C.  (r)  is  the  concentration 

is  is 

of  isotope  i  in  particle  s  at  radius  r 

f.  A  summary  of  the  amounts  of  material  in  gas  and  solid 
phases  at  the  final  temperature 


4.  3.  FUNCTION  DIFFUS  (DT,  R) 


This  function  computes  a  solution  to  the  diffusion  equation  for  a 


spherical  body  of  radius  r.  The  true  solution  is: 


,..(5) 


or 


where 


(6) 


and 


r*  i  6  V*  1  /  2  2* 

D  =  1  -  —  }  i  —  exp(-n  tt  6t/r  ) 

w  n?l  n 

D  =  6/5l\'/2L->/2+2f'ierfcJlL\.3/m 

\  r  /  \  fa  -ml  \r2) 


■r 


ierfc(x)  =  I  erfc(u)du 
'x 


erfc(x) 


,  Poo  2 
2  /  -u 

=  —  /  e 

Jx 


du  . 


(12) 


(13) 


T 

See  Appendix  D  for  sample  output. 

This  printout  and  the  calculation  are  optional  (see  Section  2). 
(See  Appendix  A  for  details.  ) 
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Function  DIFFUS  uses  the  following  approximations  to  these  solutions: 

Let  x  =  r/<s/67 


1.  If  x  >  3.  14, 


which  is  a  truncation  of  Eq.  (13)  with  the  infinite  series  neglected. 

2.  If  2.  0  <  x  <  3.  14, 

D  =  1  -  — ■  (exp(-tr*7x^)  +  —  exp(-4ir2/x2)] 
ir 

3.  If  0.  745  <  x  <  2.  0, 


D  -  1  -  ■  [exp(-ir2/x2)] 

IT 


4.  If  x  <  0.  745, 


D  =  1  , 


where  2,  3,  and  4  are  truncations  of  Eq.  (12). 

The  errors  incurred  by  these  approximations  are  bounded  as  follows: 

|E.|<  2.  5  x  10'8 
4  .5 

|e3|<  1  x  10 

|e2|<  i  x  io’5 

|E,|<  5*  IO'6 


where  E^  is  the  error  using  approximation  i. 


4.  4.  SUBROUTINE  FINAL  (AVCON,  GAS,  R,  WT,  NAME) 

Subroutine  FINAL  computes  the  concentration  of  the  last  isotope  in 
the  various  particles  assuming  that 

1.  All  other  isotopes  have  completely  decayed  into  the  last  one 

2.  None  of  the  last  isotope  has  decayed  away 


3.  All  of  the  material  in  the  gas  phase  when  the  main  process  was 
terminated  has  condensed  onto  the  particles  according  to 
surface  area 

Let  a^  =  average  concentration  of  isotope  i  in  particle  of  size  s, 

Z  =  r  /  r  where  r^  is  the  radius  of  the  largest  particles,  and 
g^  =  amount  of  isotope  i  in  the  gas  phase  when  main  process  was 
terminated. 

Then,  LC 

f 

s  /  j  is 
i=l 


G  = 


V  =  Z  w  (proportional  to  the  surface  area) 

8  8  8 

P 


where  w^  =  weight  of  particles  of  size  s. 


Also, 


Z  G 

s 

V 


+  f 


s 


and 

W  =  F  w  , 

S  8  8 

where  Fg  =  total  average  concentration  (includes  gas  condensation)  in 
particle  of  size  s 

Wg  =  total  weight  of  isotope  in  particle  of  size  s 
f  =  average  concentration  in  particle  of  size  s  due  to  the  decay  of  solids. 
The  average  concentration  of  the  last  isotope  of  the  chain,  fg,  is  considered 
to  be  equal  to  the  sum  of  the  average  concentrations  of  all  the  members  of 
the  chain  when  the  main  process  was  terminated.  Subroutine  FINAL  also 
lists  out  f,  F,  and  W  for  each  particle  size. 
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4.  5.  FUNCTION  PROFIL  (DT,  r.Q) 


Function  PROFIL  computes  the  concentration  at  a  distance  a  =  rQfrom 
the  center  of  a  spherical  particle  of  radius  r.  The  value  of  the  concentration 

/E  ) 

is  found  from  the  expression:' 


or 


.  j  +  2r^  LzAJl.  gin  £ZE£  e-n2Tr25t/r2  (14) 

C  ira  Zj  n  r 

o  n=  1 


erfcp-  +  ^  ~a]  .  erfcfeSLi  1)r  4a  1)  , 

2^/57  J  Zsfbt  I 


(15) 


where  r  =  radius  of  particle 

a  =  distance  from  center  of  particle,  a  >  0.  lr 
Cq  =  surface  concentration  (assumed  constant) 

Function  PROFIL  computes  the  value  of  the  ratio  C/Cq  using  a  truncation 
of  Eq.  (14)  or  Eq.  (15).  In  all  cases,  the  series  is  terminated  when 
either 


or  when  n  >  100,  where  Sk  is  Eq.  (14)  or  Eq.  (15)  truncated  to  k  terms. 
The  series  used  is  chosen  using  the  following  criteria: 


If  BAt  <  b.,  Eq.  (15)  is  used  where  j  is  determined  as  the 
J 

smallest  integer  such  that 

r  >  10'j, 

where  b.  =  10^  ^  )/2] 

J 

and  where  [x]  =  greatest  integer  I  such  that  I  >  x. 
Otherwise  Eq.  (14)  is  used. 


t  =  0.  01  at  present. 
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This  scheme  was  determined  by  examining  the  values  of  Eqs.  (14) 
and  (15)  for  various  numbers  of  terms  and  picking  the  equation  that 
required  the  least  number  of  terms. 


4.  6.  FUNCTION  ERF(X) 

Function  ERF  computes  an  approximate  value  of  the  error  function 
using  a  rational  approximation.  The  error  function  is  defined  as: 


erf(z)  =  —~ 
\fn 


ERF  uses  the  approximation 


(5) 


erf(x)  =  1  -  (a j  t  +  a2t 


V 


where 


and 


1  4-  px 


p  =  0.  47047 
=  0.  3480242 
a  =  -0.0958798 
a3  =  0.  7478556 


The  error  E  in  this  approximation  is  bounded  by 

|E|  <2.  5  x  10'5  . 


4.  7.  SUBROUTINE  DELTA  (N,NT,I,  PR) 

Subroutine  DELTA  computes  the  decay  coefficient,  PR,  associated 
with  decay  of  isotope  i  to  the  presently  considered  isotope  at  time  NT. 
The  array  I  specifies  the  path  being  taken,  with  1(1)  specifying  the  time  at 
which  i  is  considered. 
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A  product  of  terms  of  the  form 

A 

^inm^  "  **i,  m+l,m+l^ 

is  found.  The  d  array  is  initially  calculated  in  the  main  routine  in  such  a 
way  that 

d.  =1  if  n  >  m  . 
inm 

The  general  product  then  has  the  form 

PR  =  d.  (1  -  d.  .  ,)d  ,  { 1  -  d,  „ . .  . ,,)•••  d  .  , 

inm  i,m+l,m+l  k,m+2,£  k,£+l,  £  +  1  n,  r,  NT-1 

where  i  is  the  subscript  of  the  isotope  presently  under  consideration.  The 
values  of  n,  l,  and  r  depend  on  the  path. 


4.  8.  SUBROUTINE  BCALC  (N,  NT,  I,  KLAG,  B,  DP,  DT) 

This  routine  computes  the  correct  6t  term  associated  with  the  path 
specified  by  I  and  calls  DIFFUS  to  solve  for  a  diffusion  term. 

The  6t  term  described  by  Eq.  (A- 10)  is  computed  from  the  relation 


6t=^6  At  6..,..  At  +  .  . .  +  At  , 

Lu  iP  P  Ay  (*+l)P  P  <i+r)p  p 

p=n  p=k  p=v 

where  n.m,  k,  £,*••,  v,  and  q  are  specified  by  the  I  array.  The  corres 
ponding  diffusion  terms, 


D 


(•) 

i,  n,  m;i+l,  k,  £; 


i+r,  v,  q 


are  then  found  by  calling  DIFFUS  for  each  s.  Depending  upon  KLAG, 
either  the  first  D  or  the  second  D  in  Eq.  (A-9)  is  found  and  stored  in  the 
DD  array.  6t  is  stored  in  DT  and  the  expression 


is  stored  in  B. 


i 

s  =  1 


(D^8^  )w 

'  x,  n,  m;i+l,  k,  i+r,  v,  q'  s 


See  fection  4.  1,  item  6  for  definition  of  d 
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APPENDIX  A 


DERIVATION  OF  EQUATIONS  USED 

A.  1.  THE  GENERAL  TERM  IN  THE  SURFACE  CONCENTRATION 
EQUATION 

An  analysis  of  the  case  of  a  single  stable  isotope  with  no  precursors, 
which  eliminates  the  decay  question,  makes  it  easier  to  understand  the 
general  case.  If  Eqs.  (1)  and  (2)  in  Section  3  are  rewritten  leaving  out  the 
isotope  identification  and  letting  w^  =  1  and  p  =1,  the  following  equations 
are  obtained: 

CjGj  +  CjDj  j  =  Y  (A- 1 ) 

and 

C2G2+ClD12+(C2-Cl)D22=1f-  (A-2) 

where 

„  vnmhn 

N-  RTn  • 

Equation  (A-2)  can  be  rewritten  as  follows: 

C2C2  +  C2°22  +  C1  ®12-D22,  =  Y  ' 

The  expression  at  t^  is 

C3°3  tClD!3  +  <C2  '  C1,D23  +  lC3  '  'C2  *  Cl»  *  C,1D33  *  Y 
or 

C3G3+C1D13+<C2-C.»D23  +  (C3-C2>D33SY  ' 

The  expression  at  t^  is 

C4G4  +  Cl°14  +  (C2  -  C1)D24  +  <C3  '  C2,D34  +  <C4  '  C3)D44  =  Y  ' 
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The  expression  at  step  N  is  then 


C  G  +  C  D  + 
N  N  1  IN 


«VCk-l>DkN  *  Y  • 


(A-3) 


Equation  (A-3)  can  be  rewritten  as  follows: 


C..G,,.  +  C.Dikt  +  >  C.  D. .. 
N  N  1  IN  /  >  k  kN 


N- 1 

^  CkDk+l ,  N 


CktGki  +  >  C.  D.  Kf  -  >  ,  C,  D,  =  Y 
N  N  Lmi  k  kN  k  k+1 ,  N 


or,  finally,  as 


CN(GN  +  DNN>  +  2-Ck(DkN-Dk+l.N>=Y  ' 


(A-4) 


Next,  consider  the  general  expression  with  decay  from  precursors 
Figure  A-  1  represents  n  isotopes  at  N  time  steps;  the  connecting  lines 
represent  the  ways  in  which  decay  can  occur.  The  numbers  occurring  with 
each  state  represent  the  number  of  ways  in  which  that  isotope  at  that  time 
can  decay  to  isotope  1  at  t^.  It  is  apparent  that  if  N  <  n,  the  diagram  is 
triangular  and  is  a  Pascal's  triangle,  the  numbers  representing  binomial 
coefficients.  Thus  for  N  <  n,  the  total  number  of  terms  in  the  surface 
concentration  expression  is 

2°  +  2l  +  22  +  •  •  •  +  2N_1  =  2N  -  1  . 


For  N  >  n,  the  diagram  is  a  truncated  Pascal's  triangle.  In  this  case, 
there  are  2n  -  1  terms  from  the  purely  triangular  portion,  plus 

n-  1  TN-  1  ,  1 


n-  1  TN-1  , 

g  g(! 
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from  the  rectangular  portion,  giving  a  total  of 


n- 1 


+  2 


n 


1  . 


( A-  5 ) 


In  Eq.  (A-4),  each  term  in  the  sum  consists  of  a  surface  concentration 
times  a  diffusion  term.  In  each  case  the  diffusion  term  is  the  difference 
of  two  diffusion  terms,  one  representing  diffusion  from  t^  to  t^  and 
the  other  from  t^+jj  t0  ^  Each  of  the  M  terms  in  the  general  expression 
is  of  a  similar  form,  modified  by  a  decay  term. 

As  an  example,  consider  the  terms  associated  with  the  position  n  =  2 
in  Fig.  A.  1  at  t^.  There  are  three  possible  paths  for  decay  to  isotope  1. 

*  The  decay  coefficients  associated  with  these  paths  are 


(1  '  d21 1  >dI23  ] 
d211(1  ‘  d222>d133/ 
d212(I  *  d233>  ) 


where  d^^represents  the  decay  coefficient  of  isotope  i  from  time  step 
Atn  through  At^.  The  corresponding  isotope  indexed  diffusion  terms  are: 

°2,  1,  1;  1 , 2,  4  "  Dl,2,4 
D2,  1,  2;1,  3,  4  ‘  °2,  2,  2;  1 ,  3,4 
°2,  1,  3;  l ,  4,  4  ’  °2,2,  3;  1 , 4,  4 


where  D.  .  ,  ,  .  represents  diffusion  as  isotope  i  from  t  to 

i,  n, m;  l+l, k,£; . . .  ;  i+r ,  v,  q  r  r  n 

t  »  as  isotope  i+1  from  t,  tot.,  etc.,  and  as  isotope  i+r  from  time  t  to  t  ^ 
m  r  k  t,  r  v  q 

(see  Eq.  (A-10)). 


See  Section  4.  1,  item  6. 

TNote  that  decay  is  carried  only  through  Step  3,  while  diffusion  is 
carried  through  Step  4  (see  Section  3). 


The  total  contribution  to  isotope  1  at  t^  from  isotope  2  at  tj  can  be 
written  as 


'21 


'  d2 1 1  )dl 2 3<D2,  1.  1 ;  1 , 2,  4 


1,2,4 


) 


d211(I  ‘  d222)d133(D2,  1.2;1,  3,4 


°2,  2,  2;  1 ,  3,  4* 


+  d212(1  '  d233,(D2,  1,  3;  1,4,  4  "  °2,2,  3;  1,4,  4^ 


(A-8) 


where  again  a  single  particle  of  unit  weight  is  assumed. 

From  the  above  discussion  it  can  be  seen  that  the  general  term  is 
of  the  form 


Cindinm^  ‘  di,m+l.m+Pdi+l,m+2. k^1  "  d  i+1,  k+1 ,  k+1  ^ 


d  i+r,  l,  N-  1 


(D 


(a) 

i.n,m+l;i+l.m+2.k+l 


i+r,  l,  N 


where 


-D 


(a) 

i,n+l,m+l;i+l,m+2,  k+1; 


i+r,  t,,  N 


)w. 


; 


(A-9) 


1.  We  are  considering  the  contribution  to  isotope  i  +  r  at  time  N 
via  some  path,  from  isotope  i  at  time  n. 

2.  m  >  n  +  1 

k  >  m  +  1 


l  >  i  +  r 

3.  There  are  p  particle  sizes  with  weights  Wj,  W£>  ....  Wp. 

4.  If  n  >  m,  djnm  =  1  as  a  mathematical  convenience. 

The  quantity  d!8*  . -  represents  a  solution  of  the 

ii  Hi  rtif l *  1 1  K|  if i  i  iiTf  Vf  cj 

diffusion  equation  for  a  particle  of  radius  rg  and  6t  term  given  by 

n  i  q 

5t=Y^5.  At  +  V*  5...  At  +  •  •  ■  +V5U  At  (A- 10) 

in  n  i+I,n  n  l+r,n  n  '  ' 

n=q  n=K  n=v 
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where  6^  is  defined  by  Eq.  (11),  Section  4.  1,  item  8,  and  Atn  is  defined  by 
Eq.  7,  Section  4.1,  item  4. 


A.  2.  AVERAGE  CONCENTRATIONS  AND  PROFILES 


As  in  the  case  of  surface  concentrations,  the  case  considered  is 
that  of  a  single  stable  isotope  with  no  precursors.  The  basic  relationship 


(A-  1 1 ) 


where  C  is  an  average  concentration  increment.  C  is  a  surface  concen¬ 
tration  increment,  and  Disa  diffusion  term  as  in  Section  A.  1.  At  the 


first  time  step,  it  is  clear  that 


ci=ciDn  ' 


At  the  second  time  step, 


c2  =  CIDia  +  <c2  -  C1>D22  ' 

where  the  first  term  on  the  right  represents  the  concentration  at  time 
step  1  further  diffused  and  the  second  term  represents  material  condensed 
and  diffused  during  time  step  2.  Similarly,  at  time  step  3, 

^3  =  CIDI3  +  (C2  *  C1>D23  +  <C3  -  C2»D33  ■ 
where  the  three  terms  represent,  respectively,  condensation  during  time 
steps  1,  2,  and  3,  and  diffusion  from  those  times  to  time  3.  At  time 


step  N,  it  follows  that 


CN  C1D1N  + 


(Ck  '  Ck-l)Dk,  N  * 


(A- 12) 


where  each  term  in  the  sum  represents  a  concentration  diffused  over  the 
specified  time  interval. 


34 


s 


To  obtain  an  expression  similar  to  Eq.  (A-4),  Eq.  (A-12)  is  rewritten 


ab 


N  - 1 


CKI  =  C,D1KT  +  (CKI  -  C„  ,)Dkt  „t  +  23  C.  d 


N-l 


N  ''l'  lN  ’  '~N  ''N-l'^N.N  ’  ~k  k,  N 

k=2  m 


Ck-lDk,  N  ’ 


or 


or,  finally,  as 


N-l  N-l 

'N  =  CNDN,N+E  CkDk,  N  ‘  ^3  CkDk+l,N 
N-l 

CN  =  CNW£  Ck^Dk,  N  ■  Dk+l,N)  ‘ 


(A-  13) 


It  can  be  seen  that  the  general  expression  including  decay  from 
precursors  consists  of  terms  similar  to  Eq.  (A-9).  The  number  of  such 
terms  is  the  same  as  that  described  in  Section  A*l,  since  there  is  a  con¬ 
tribution  from  each  isotope  at  each  preceeding  time  step  for  each  possible 
decay  path.  Here  the  expression  for  average  concentration  in  the  s**1 
particle  size  consists  of  terms  of  the  form: 

Cindinm^  "  di,m+l,m+P  di+l,  m+2,k^  ’  d i  +  1 ,  k+1 ,  k+1 ^ 

(s) 

d  i+r,  l,  N-  1  ^i,n»m+l;i+l ,  n+2,  k+1;  •••  ;  i+r,£,N 

-  d(8) 

i,n+l,m+l;i+l ,  m+2,  k+1;  •  •  •  ;  i+r,  l,  N)  (A-14) 
summed  over  all  contributions. 

To  find  profiles,  it  is  desired  to  find  the  concentration  at  specific 
points  within  the  particle.  The  basic  relation  here  is  (see  Ref.  6) 


=  F(a,  r,  6t)  , 


(A-  15) 


,  Tt  ,  ,  2r  (-1  )n  .  mra  ,  2  2.  ,  2 

where  r(a,  r,  6t)  =  1  +  —  )  M  sin -  exp(-  n  tr  6t/r  )  . 

na  H  nr 
n=  1 
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Aa  with  average  concentration*,  Eq.  (A- 15)  gives  a  concentration  increment 
due  to  diffusion  during  the  given  time  and  with  zero  initial  concentration. 

For  a  stable  isotope  with  no  precursors  we  can,  in  an  analagous  manner, 
obtain  the  expression 

CiN*a)  =  CiNr(a,r'6iNAtN)  Cik  [r(a’  r’  ^6inAtn)  ’  ^  r’  ^  6inAtn)] 

n=  +1  (A16) 

and  in  the  general  case  obtain  a  sum  of  terms  similar  to  Eq.  (A.  14)  with 
each  D  replaced  by  the  appropriate  T. 
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APPENDIX  B 


MODEL  CALCULATIONS 


Some  of  the  initial  results  of  calculations  obtained  using  this  program 

are  considered  in  this  appendix.  In  these  calculations  examples  of  a  refractory 

and  a  volatile  chain,  95  and  137,  respectively,  were  selected  for  study.  The 

results  should  demonstrate  the  difference  in  behavior  for  extreme  cases  of 

refractory  and  volatile  nuclide  chains  with  varying  conditions  and  thus  pro* 

vide  limits  of  fission  product  behavior  during  fallout  formation.  Data  sources 

for  these  calculations  have  been  listed  earlier.  In  addition  to  these  sources, 

the  standard  particle  distribution  shown  in  Table  B.  1  was  developed  from 

* 

one  used  by  Miller,  experimental  and  estimated  diffusion  coefficients  were 
employed  as  measured  or  estimated  in  reported  studies,  and  total  fissions 
were  scaled  as  a  function  of  weapon  fission  yield  according  to  Miller.  ^ 

Several  program  variables  had  to  be  investigated  in  order  to  demon¬ 
strate  the  degree  of  dependence  of  the  calculational  system  on  these  variables. 
For  instance,  the  variation  of  results  with  the  size  of  the  temperature  incre¬ 
ment  must  be  vanishingly  small  for  the  selected  temperature  increment  as 
smaller  increments  are  considered.  Also,  the  variation  of  results  with 
choice  of  initial  and  final  calculational  temperatures  should  be  defined. 

Finally,  the  fact  that  the  particle  size  distribution  employed  is  idealized— 
that  is.  some  selected  radii  represent  all  of  the  mass  of  the  particles— should 
be  considered.  These  points  have  been  investigated  and  seem  to  be  satisfac¬ 
torily  understood  as  will  be  demonstrated  as  the  calculations  are  presented. 

137  1 3  () 

The  calculations  for  the  95  and  137  chains  (neglecting  I  -*  Xe  and 

1*^®  -*  Xe*^)  were  performed  using  2700°K  as  an  initial  temperature.  200°K 

increments,  terminating  at  900°K,  and  employing  the  particle  distribution 

described  in  Table  B.  1.  The  size  of  the  detonation  was  taken  as  100  KT 

9 

pure  fission.  The  amount  of  soil  associated  with  the  cloud  was  8.  5  *  10  g. 

$ 

Miller.  C.  F.  ,  private  communication,  November  1965. 
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Table  B.  1 

STANDARD  PARTICLE  SIZE  DISTRIBUTION 


Designation 

Particle  Radius 
(cm) 

Weight  Percent  of 
Particle  Fraction 

T 

0.  00030 

0.  38 

S 

0.  0030 

28.  92 

M 

0.  0154 

60.  2 

L 

0.  13 

10.  5 

The  calculated  descriptions  of  behavior  of  the  137  and  95  chains  are 
presented  in  Figs.  B.  1  and  B.2,  respectively.  In  these  figures,  theaverage 
concentration  of  a  particular  nuclide  in  the  chain  is  presented  as  a  function 
of  the  last  temperature  for  which  the  ad  sorption- diffusion  process  was  cal¬ 
culated.  Also  included  are  the  gas  phase  contents  of  the  nuclides  at  the 
same  temperatures  (designated  by  open  circles).  In  the  phenomena  des¬ 
cribed  in  Fig.  B.  1,  the  curve  for  the  tellurium  gas  phase  content  is 
decreasing  as  the  result  of  two  processes: 

1.  Te^37  decay. 

2.  Absorption  of  Te^7  jn  the  fallout  particles. 

Corresponding  processes  are  involved  to  different  degrees  in  the  description 
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of  the  gas  phase  content  of  any  nuclide.  In  the  case  of  Te  at  the  lowest 

137 

temperature,  so  little  of  the  Te  has  been  absorbed  that  the  gas  phase 
content  curve  closely  represents  the  tellurium  decay.  The  gas  phase  con¬ 
tent  curves  for  iodine  and  xenon  also  essentially  represent  radioactive 
decay  phenomena  alone,  since  both  of  these  elements  are  very  volatile. 

On  the  other  hand,  the  gas  phase  content  of  cesium  is  dropping  sharply  with 

137 

temperature,  not  because  Cs  is  decaying  but,  in  this  case,  because  94% 

of  the  cesium  in  the  system  has  been  absorbed  by  the  fallout  particles  by 
o 

the  1900  K  point  and  more  is  absorbed  at  lower  temperatures.  Thus,  the 
cesium  curve  describes  mainly  the  absorption  process. 

The  average  concentrations  of  Te^^  as  they  vary  with  temperature 
and  particle  size  were  also  investigated.  The  large  particles,  L,  as 
defined  in  Table  B.  1,  shew  a  slow  but  continuous  increase  in  average 

tellurium  concentration  as  the  temperature  decreases  (see  Fig.  B.  1). 
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AVERAGE  CONCENTRATION  (G/G) 


Fig.  B.  [--Calculated  absorption  in  fallout  of  members  of  137  chain 
in  g  fission  product  per  g  fallout  as  a  function  of  temperature  and 

particle  size 
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GAS  PHASE  CONTENT  (G) 


AVERAGE  CONCENTRATION  (G/G) 


Fig.  B.  2- -Calculated  absorption  ir  fallout  of  members  of  95  chain 
in  g  fission  product  per  g  fallout  as  a  function  of  temperature  and 

particle  size 
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GAS  PHASE  CONTENT  (G) 


Also,  at  the  initial  temperature  the  L  particles  exhibit  a  slightly  lower 
concentration  than  the  medium* sized  particles,  M.  This  point  is  important 
because  it  suggests  that  at  the  highest  temperature  of  the  calculation,  the 
L  particles  did  not  equilibrate  with  the  gas  phase;  that  is,  diffusion  was 
not  complete.  In  contrast,  the  tiny  particles,  T,  the  small  particles,  S, 
and  the  M  particles  exhibit  essentially  the  same  average  concentration 
at  the  initial  temperature.  Further  consideration  of  Fig.  B.  1  reveals  that 
the  L  particles  experience  a  somewhat  smaller  increase  in  tellurium 
concentration  than  the  M  particles  as  the  temperature  decreases.  Also, 
the  S  and  M  curves  exhibit  a  departure  from  the  more  sharply  rising  T 
curve  at  two  different  lower  temperatures.  Departure  points  from  the  T 
curve  occur  as  the  central  portions  of  the  L,  S,  a  i  M  particles  effectively 
lose  contact  with  the  surface  and  diffusive  transport  begins  to  control  the 
rate  of  absorption.  In  addition,  at  1300°K  absorption  into  the  T  particles 
is  diffusion-controlled  (at  this  temperature  the  surface  concentration  is 
calculated  to  be  twice  the  average  concentration). 

The  behavior  of  other  elements  in  the  137  chain  can  also  be  observed 
in  Fig.  B.l.  Neither  iodine  nor  xenon  is  strongly  absorbed  at  even  the 
lowest  temperature  of  this  study.  The  higher  temperature  portions  of  the 
iodine  and  xenon  curves  in  Fig.  B.l  reflect  the  low  absorption  and  a  low 
temperature  coefficient  of  absorption  for  these  isotopes.  At  lower  tem¬ 
peratures,  however,  the  iodine  curves  reflect  the  behavior  of  the  parent 
tellurium  curves.  That  is,  at  these  temperatures  the  main  iodine  concen¬ 
tration  in  the  particles  comes  from  absorbed  tellurium  decay  to  iodine,  whose 
loss  from  the  particle  is  diffusion- controlled.  The  low  temperature 
increase  in  xenon  concentration  in  the  particles  is  due  to  iodine  decay, 
and  xenon  loss  is  also  diffusion-controlled.  This  latter  process  has  just 
become  apparent  as  a  result  of  the  lowest  temperature  calculations  shown 
in  Fig.  B.  1.  Similar  phenomena  can  be  observed  in  the  cesium  absorption 
curves.  Fission  cesium  is  mainly  absorbed  at  2100°K,  where  the 
absorption  curves  in  Fig.  B.l  stop  their  steep  ascent.  This  is  the  same 
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temperature  region  where  gas  phase  cesium  starts  its  sharp  decrease. 

Note  that  the  concentration  of  cesium  in  the  L  particles  does  not  appre¬ 
ciably  change  at  lower  temperatures.  The  remaining  gaseous  cesium,  and 
that  formed  in  the  gas  phase  decay  of  xenon,  is  picked  up  mainly  by  the  smaller 
particles  because  of  their  large  surface  to  volume  ratio.  The  cesium 
concentration  of  the  M  particles  does  not  change  much  after  1900°K, 
while  that  of  the  S  particles  begins  to  level  off  at  1700°K.  The  cesium 
concentration  of  the  T  particles  increases  sharply,  reflecting  the  continued 
absorption  of  cesium  produced  from  xenon  (g). 

l'he  behavior  of  the  137  nuclide  chain  can  thus  best  be  explained  in 
terms  of  diffusion-controlled  sorption  of  tellurium  and  cesium  and  decay 
of  the  nuclides.  Note  that  the  sorption  processes  depicted  in  Fig.  B.  1 
account  for  only  a  minor  portion  of  the  total  mass  of  fission  products  of 
this  nuclide  chain.  The  rest  of  the  material,  at  times  subsequent  to  tho^c 
of  this  calculation,  will  go  through  the  process  of  decaying  to  cesium  and 
depositing  on  the  outer  portions  of  the  fallout  particles  or  on  any  condensed 
phase  with  which  it  comes  in  contact.  In  this  study,  however,  it  is  assumed 
that  the  material  remaining  in  the  gas  phase  at  900°K  is  surface-deposited 
on  the  fallout  particles  themselves,  according  to  their  surface  area.  In 
reality,  xenon,  with  a  relatively  long  half  life  (234  sec),  may  well  separate 
from  portions  of  the  fallout  particles  before  decaying  and  depositing. 

Figure  B.2  shows  similar  data  for  the  95  chain.  The  behavior 
depicted  in  this  figure,  except  for  that  of  rubidium, represents  essentially 
complete  absorption  at  the  highest  temperature  of  the  study,  2700°K,  and 
rubidium  absorption  is  essentially  complete  by  1900°K.  The  absorption 
differences  between  the  smaller  particles  and  the  L  particles  result  from 
the  nuclides  in  this  chain  being  incompletely  diffused  in  the  larger  particles 
at  2700°K.  The  splitting  if  the  M,  S,  and  T  curves  of  rubidium  results 
from  diffusion  control  becoming  important  at  lower  temperatures  for 
smaller  particles,  as  was  discussed  previously  for  tellurium  in  the  137 
chain. 
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A  further  description  of  the  effect  of  initial  temperature  is 
appropriate  here.  Table  B. 2  contains  the  900°K  average  concentration  ratios 
of  L  and  M  for  zirconium,  yttrium,  and  strontium  and  the  L/M,  M/S, 
and  S/T  average  concentration  ratios  for  rubidium.  For  a  highly  absorbed 
species,  the  ratio  L/M  differs  from  unity  but  obviously  will  approach  unity 
as  starting  temperatures  are  increased.  However,  for  rubidium,  which 
is  not  so  highly  absorbed,  the  concentration  ratios  in  various  sized 
particles  are  independent  of  the  starting  temperature;  and  the  smaller 
the  particles,  the  more  nearly  this  value  approaches  unity.  The  cal- 
culational  values  for  less  highly  absorbed  species  thus  appear  to  be 
independent  of  starting  temperature. 


Table  B.  2 

CONCENTRATION  RATIOS  AT  900°K  OF  PARTICLES  OF  DIFFERENT 
SIZE  FOR  VARIOUS  STARTING  TEMPERATURES 


Starting 

Temperature 

(°K) 

Average  Concentration  Ratio 

Zr  (L/M) 

Y  (L/M) 

Sr  (L/M) 

Rb  (L/M) 

Rb  (M/S) 

Rb  (S/T 

2700 

0.  378 

0.  497 

0.61 1 

0.  398 

0.  844 

0.  950 

3100 

0.600 

0.  754 

0.  859 

0.  398 

0.  844 

0.  950 

3900 

0.  862 

0.  924 

0.  883 

0.  398 

0.  844 

0.  950 

The  phenomenon  of  lack  of  equilibration  of  nonfractionated  nuclides 
in  the  larger  particles  at  high  temperatures  which  is  seen  when  the  present 
model  is  used  would  seem  to  indicate  some  fractionation  with  size.  An 
investigation  of  fallout  should  reveal  whether  this  is  indeed  the  case. 

The  orientation  calculation  described  above  has  demonstrated  many 
of  the  features  of  the  diffusion-controlled  model.  These  features  seem  to 
find  some  support  in  studies  of  fallout  and  certainly  should  be  studied  more 
thoroughly.  However,  before  further  questions  concerning  calculations 
of  fallout  phenomena  are  investigated,  some  of  the  strictly  mathematical 
variations  within  the  model  calculations  should  be  considered. 
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An  inveitigation  of  different- sized  temperature  stepping  intervals 
was  made.  In  general,  in  the  calculations  reported  here,  200°K  tem¬ 
perature  intervals  were  used.  However,  these  test  calculations,  which 

o  o 

were  initiated  at  2700  K  and  were  carried  through  900  K,  were  made  for 

the  two  nuclide  chains  in  both  200°  and  100°K  increments.  Comparing 

results  was  more  complicated  than  it  first  appeared  because  the  final 

o  o 

diffusion  time  corresponded  to  800  K  in  the  200  interval  calculations  and 

to  850°K  in  the  100°  interval  calculations.  In  order  to  avoid  this  difficulty, 

both  systems  were  allowed  to  decay  and  to  absorb  the  resulting  gaseous 

1  37 

fission  products  (mainly  Cs  )  according  to  surface  area,  thereby  nearly 
eliminating  the  final  calculational  temperature  differences.  The  com¬ 
parison  of  results  is  presented  in  Table  B.3.  The  only  difference  that  seems 
significant  occurred  in  the  L  particle  zirconium  concentrations,  and  it 
can  apparently  be  explained  by  the  difference  in  initial  starting  temperatures. 
The  particles  had  essentially  twice  the  time  to  diffuse  at  2700°K  in  the  case 
of  the  200°  increments  as  in  the  case  of  the  100°  increments.  While  a  better 
example  might  be  constructed,  it  would  appear  from  this  study  that  the 
size  of  the  temperature  increment  used  in  the  model  may  not  be  very 
important  if  chosen  as  200°K  or  less. 


Table  B.  3 

EFFECT  OF  SIZE  OF  TEMPERATURE  INCREMENT  ON 
FISSION  PRODUCT  ABSORPTION 


Nuclide 

Particle 

Size 

Absorption  (g) 

100°C  Temp. 
Interval 

200°C  (Std) 

Temp.  Interval 

Zr 

T 

0.456 

0.  456 

S 

34.67 

34.  63 

M 

71. 52 

71.  17 

L 

6.  86 

7.  26 

Cs 

T 

14.62 

14.63 

112.  40 

1 12.  44 

1 

46.  92 

46  88 

1.045 

1.  042 

44 


r 


In  the  investigation  of  the  effect  of  the  particle  size  distribution,  three 
different  distributions  were  studied  under  the  same  conditions.  The  distri¬ 
butions  employed  and  the  weights  of  soil  associated  with  each  distribution  in 
this  case  are  presented  in  Table  B.4.  The  amounts  of  cesium  and  zirconium 
fission  products  (in  grams)  after  decay  of  precursor  fission  products  are  also 
presented  for  each  fallout  particle  size  as  in  the  previous  case.  Finally, 
a  ratio  of  the  cesium  and  zirconium  average  concentrations  associated 
with  each  particle  size  is  presented.  This  ratio  is  proportional  to  a 
fractionation  factor  ratio  comparing  these  elements.  The  proportionality 
constant  is  0.  649,  a  nuclide  chain  yield  ratio.  Further  consideration  of 
Table  B.  4  reveals  that  fractionation  factors  for  the  very  smallest  particles 
far  exceed  unity,  while  for  the  larger  particles  these  factors  are  less  than 
unity.  Fractionation  factors  in  this  model  are  thus  quantities  that  depend 
on  particle  size  (and  also  probably  on  the  degree  of  detachment  of  the  fallout 
particles  from  the  nuclear  cloud  before  conversion  of  volatile  fission 
products  to  condensable  ones).  Therefore,  if  this  model  is  correct, 
observed  fractionation  factors  should  be  dependent  on  particle  size  and 
thus  not  accurately  definable  without  ample  description  of  this  parameter. 
Indeed,  these  factors  should  be  sensitive  not  only  to  the  particle  size  in 
question  but  also  to  the  particle  size  distribution.  Specifically,  one 
observes  using  distributions  (2)  and  (3)  that  fractionation  factors  for  the 
same  size  particles  can  differ  considerably.  Particle  distribution 
(3)  contains  a  large  amount  of  very  small  particles.  Distributions 

* 

(1)  and  (2)  are  similar,  since  they  both  are  based  on  Miller's  description. 

Note  that  the  distributions  (1)  and  (2)  produce  similar  fractionation  factors 
for  similarly  sized  particles.  It  seems  apparent  that  in  testing  this  model 
for  reproduction  of  fractionation  factors,  one  must  specify  the  particle 
size  distribution,  a  factor  often  not  well  documented. 

A  study  was  made  of  the  effect  of  the  size  of  the  detonation  and  the 

I 

amount  of  soil  in  the  cloud  (that  is,  the  proximity  of  the  detonation  to  the 
* 

Miller,  C.  F. ,  private  communication,  November  1965. 
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land  surface).  The  calculated  final  amounts  of  cesium  and  zirconium  for 
the  various  particle  sizes  as  a  function  of  these  two  quantities  are  shown 
in  Table  B.  5.  Note  that  the  calculated  fission  product  distributions  in  the 
particles  were  only  slightly  different  for  the  various  conditions.  Specific 
activities  per  gram  of  fallout  vary  directly  with  yield  and  indirectly  with 
the  amount  of  soil,  but  the  relative  fission  product  distribution  among  the 
different-sized  particles  is  not  greatly  affected  by  these  parameters. 

This  reemphasizes  the  importance  of  particle  size  distributions  and 
chemistry  in  determining  the  fate  of  fission  products  during  fallc  * 
fo  rmation. 


Table  B.  5 


VARIATION  OF  FISSION  PRODUCT  ABSORPTION  WITH  SIZE 
AND  HEIGHT  OF  DETONATION  EFFECTS  FOR  A 
1 00-  KT  FISSION  YIELD  DETONATION 


Cesium  Absorbed  (g) 

Zirconium  Absorbed  (g) 

Particle 

|  | 

10  MT. 

10MT, 

Radius 

Total 

Total 

(cm) 

Standard 

1/100  Soil 

Yield 

Standard 

1/100  Soil 

Yield 

0.  0003 

14.  63 

14.  76 

13.  84 

0.  444 

0.  539 

0.432 

0.  003 

112.  44 

113.20 

1 10.  21 

33.  73 

35.  54 

32.  86 

0.  0154 

46.  88 

46.  06 

49.  55 

69.  28 

67.  76 

68.  41 

0.  130 

1.  042 

0.  954 

1.  387 

10.  06 

9.67 

11.  81 

The  ce">nuter  program  can  also  provide  fission  product  profiles. 
Calculated  profiles  for  both  the  137  and  95  chains  have  been  investigated 
for  a  100-KT  fission  detonation,  using  200°  temperature  increments  and 
the  particle  size  distribution  given  in  Table  B.  1.  The  calculations  were 
iritiated  at  ?700°K  for  the  137  chain  and  at  3100°K  for  the  95  chain.  The 
900°K  profiles  for  the  95  chain  are  presented  as  a  function  of  particle  size 
in  Fig.  B.  3.  and  those  for  the  137  chain  are  presented  in  Figs.  B.  4.  B.  5.  and 
B.6.  An  unusual  abscissa  scale  is  used  in  these  figures  solely  for  the 
purpose  of  demonstrating  the  calculated  values  more  efficiently. 
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CONCENTRATION  (G/G) 


Fig.  B. 3-  -Concentration  profiles  at  900°K  of  95  nuclide  chain 
members  as  a  function  of  particle  size  for  standard  detonation 
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Fig.  D.  4--Tcllurium - 1 37  concentration  profiles  at 
900°K  for  standard  particln  size  distribution  and 
standard  detonation 


0  10°  10  lO--*  10“*  5*10“*  10“A  0.3  0.6  0.9 

PENETRATION  (PARTICLE  RADII) 

Fig.  B.  5- -Iodine  and  xenon- 137  concentration  profiles  at  900°K 
for  standard  particle  size  distribution  and  standard  detonation 


CONCENTRATION  (G/G) 


Fig.  B.  6- - Cesium- 1 37  concentration  profiles  at  900°K  for 
standard  particle  size  and  standard  detonation 
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The  900°K  profiles  for  the  95  chain  are  very  simple.  They  show  only 

small  deviations  from  a  flat  profile  for  rubidium  in  the  S.  M,  and  L 

particles  and  for  the  other  fission  products  just  in  the  L  particles.  The 
o 

900  K  profiles  for  the  137  chain  are  structurally  more  interesting.  The 
distribution  of  tellurium  within  the  standard  particles  at  this  temperature 
is  shown  in  Fig.  B.4.  Calculated  tellurium  concentrations  span  eight 
orders  of  magnitude  in  the  L  particles  and  five  orders  of  magnitude  in  the 
T  particles.  The  concentrations  have  been  drawn  as  a  continuous  curve 
in  these  diagrams  as  an  aid  to  the  reader  and  arc  not  to  be  used  in  inter¬ 
polating.  The  tellurium  profiles,  then,  are  characteristic  of  a  moderately 
volatile  fission  product  (1%  condensed  at  900°K  in  this  calculation)  with  a 
moderately  high  Henry's  law  temperature  coefficient. 

The  iodine  and  xenon  profiles  (shown  in  Fig.  B.  5)  are  even  more 
interesting,  since  they  exhibit  maxima  and  minima  as  a  function  of  radius 
and  particle  size.  The  sharp  maxima  in  the  iodine  curves  can  be  attributed 
to  absorbed  tellurium  decay  and  subsequent  diffusion  and  vaporization  of  the 
resulting  iodine.  The  maxima  in  the  xenon  curves  can  be  attributed  to 
absorbed  iodine  decay  and  diffusion  and  vaporization  of  resultant  xenon. 

The  cesium  profiles  are  shown  in  Fig.  B.  6.  This  set  of  curves  would 
appear  to  be  exhibiting  fine  structure  that  is  somewhat  dependent  on  the 
calculational  model.  The  near- surface  behavior  of  calculated  profiles 
would  seem  to  have  an  extraneous  structure  that  is  a  function  of  the  size 
of  the  time -temperature  steps.  The  use  of  shorter  time  steps  should 
smooth  out  this  profile  considerably.  However,  the  general  structure  of 
these  profiles  seems  reasonable. 

It  is  interesting  to  consider  the  penetration  depths  used  in  these  mathe- 

-5 

matical  profile  calculations.  The  particle  radius  0. 0003  cm  times  10  fractional 

o 

penetration  gives  0.  3  A,  a  value  too  low  on  the  atomic -size  scale  to  be  a 

credible  penetration  depth  This  can  also  be  said  for  the  0.  0003-cm 

-4 

radius  particles  at  10  radius  penetration,  and  the  0.  003-cm  radius 
particles  at  10  ^  radius  penetration.  A  value  between  15  and  30  X,  a 
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borderline  penetration,  ia  the  next  smallest  penetration  depth  plotted. 

The  subsequent  minimum  value  is  130  X.  a  value  such  that  diffusion  is 
reasonable  to  consider. 

A  plot  of  the  iodine  profiles  for  the  0.  0003-cm  radius  particles,  T, 
as  a  function  of  temperature  is  shown  in  Fig.  B.  7.  These  curves  demon* 
strate  a  general  decrease  in  total  iodine  loading  as  the  cloud  cools.  This 
can  be  attributed  to  the  expansion  of  the  cloud,  high  iodine  diffusion 
coefficients,  and  a  low  iodine  Henry's  law  temperature  coefficient.  At 
the  lower  temperatures,  where  iodine  is  not  very  mobile,  the  iodine 
profile  reflects  the  parent  tellurium  profile.  The  difference  between  the 
1100°K  and  the  900°K  profiles  ie  caused  by  additional  tellurium  decay 
to  iodine. 

This  orientation  investigation  demonstrates  the  effects  of  condensed 

state  diffusion  in  a  fallout  calculational  model  embodying  this  phenomenon. 

It  is  interesting  to  compare  this  system  with  the  Miller  model,  where 

o 

diffusion  is  deemed  infinitely  fast  above  1673  K  and  negligibly  slow  below 
this  temperature.  For  the  case  of  the  highly  absorbed  nuclide  chain,  the 
only  noticeable  difference  is  that  the  calculations  for  the  larger  particles 
in  the  diffusion-controlled  model  suggest  that  equilibrium  conditions  may 
not  be  attained  by  these  particles  during  fallout  formation.  Equilibrium 
conditions  are  assumed  by  Miller.  In  the  case  of  a  volatile  chain,  the 
diffusion  model  predicts  an  absorption  behavior  which  varies  with  particle 
size.  Miller's  model  assumes  no  absorption  difference  with  particle  size. 

In  the  diffusion  model,  cooling  rate  will  slightly  affect  absorption,  while 
in  Miller's  model  it  will  affect  absorption  only  to  the  extent  that  this  quantity 
interacts  with  fission  product  decay.  Radial  concentration  gradients  are 
predicted  within  a  fallout  particle  with  the  diffusion  model,  while  with 
Miller's  model  there  are  no  concentration  gradients  within  a  particle, 
only  at  the  surface  of  the  particle.  It  seems  reasonable  to  state  that  the 
diffusion  model  provides  a  more  revealing  and  perhaps  a  more  accurate  cal¬ 
culational  tool  for  considering  fallout.  It  is,  however,  a  limited  model,  and 
the  assumptions  underlying  the  development  of  the  model  should  not  be  ignored. 
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APPENDIX  D 


SAMPLE  OUTPUT  AND  SUBPROGRAM  LISTINGS 
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